home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / histogram / test2d_resample.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-11-01  |  2.9 KB  |  118 lines
  1. /* histogram/test2d_resample.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <math.h>
  23. #include <gsl/gsl_histogram2d.h>
  24. #include <gsl/gsl_test.h>
  25. #include <gsl/gsl_ieee_utils.h>
  26.  
  27. #include "urand.c"
  28.  
  29. int
  30. main (void)
  31. {
  32.   size_t i, j;
  33.   int status = 0;
  34.   double total = 0;
  35.   size_t N = 200000;
  36.  
  37.   gsl_histogram2d *h;
  38.  
  39.   gsl_ieee_env_setup ();
  40.  
  41.   h = gsl_histogram2d_calloc_uniform (10, 10, 0.0, 1.0, 0.0, 1.0);
  42.  
  43.   for (i = 0; i < 10; i++)
  44.     {
  45.       for (j = 0; j < 10; j++)
  46.     {
  47.       double w = 10.0 * i + j;
  48.       total += w;
  49.       gsl_histogram2d_accumulate (h, 0.1 * i, 0.1 * i, w);
  50.     }
  51.     }
  52.  
  53.   {
  54.     gsl_histogram2d_pdf *p = gsl_histogram2d_pdf_alloc (10,10);
  55.  
  56.     gsl_histogram2d *hh = gsl_histogram2d_calloc_uniform (20, 20,
  57.                               0.0, 1.0,
  58.                               0.0, 1.0);
  59.  
  60.     gsl_histogram2d_pdf_init (p, h);
  61.  
  62.     for (i = 0; i < N; i++)
  63.       {
  64.     double u = urand();
  65.     double v = urand();
  66.     double x, y;
  67.     status = gsl_histogram2d_pdf_sample (p, u, v, &x, &y);
  68.     status = gsl_histogram2d_increment (hh, x, y);
  69.       }
  70.  
  71.     status = 0;
  72.     for (i = 0; i < 20; i++)
  73.       {
  74.     for (j = 0; j < 20; j++)
  75.       {
  76.         double z = 4 * total * gsl_histogram2d_get (hh, i, j) / (double) N;
  77.         size_t k1, k2;
  78.         double ya;
  79.         double x, xmax, y, ymax;
  80.  
  81.         gsl_histogram2d_get_xrange (hh, i, &x, &xmax);
  82.         gsl_histogram2d_get_yrange (hh, j, &y, &ymax);
  83.  
  84.         gsl_histogram2d_find (h, x, y, &k1, &k2);
  85.         ya = gsl_histogram2d_get (h, k1, k2);
  86.  
  87.         if (ya == 0)
  88.           {
  89.         if (z != 0)
  90.           {
  91.             status = 1;
  92.             printf ("(%d,%d): %g vs %g\n", (int)i, (int)j, z, ya);
  93.           }
  94.           }
  95.         else
  96.           {
  97.         double err = 1 / sqrt (gsl_histogram2d_get (hh, i, j));
  98.         double sigma = fabs ((z - ya) / (ya * err));
  99.         if (sigma > 3)
  100.           {
  101.             status = 1;
  102.             printf ("%g vs %g err=%g sigma=%g\n", z, ya, err, sigma);
  103.           }
  104.           }
  105.       }
  106.       }
  107.  
  108.     gsl_histogram2d_pdf_free (p) ;
  109.     gsl_histogram2d_free (hh) ;
  110.     
  111.     gsl_test (status, "gsl_histogram2d_pdf_sample within statistical errors");
  112.   }
  113.  
  114.   gsl_histogram2d_free (h) ;
  115.  
  116.   exit (gsl_test_summary ());
  117. }
  118.